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The Pierre Auger Observatory, recently completed, has been operational since 2004. As a hybrid experiment, it 
allows for a wide range of measurements of UHECR-induced extensive air showers (EAS), including measurements 
of the EAS particle content on ground which is sensitive to high-energy hadronic interactions. We present the 
results of several independent measurements of the EAS muon content on ground in Auger data at a primary 
energy of 10 EeV. We discuss implications on high-energy hadronic interaction models and cosmic ray composition. 



1. Introduction 

Cosmic ray interactions at ultra-high energies 
offer unique insights into particle collisions at 
center-of-mass energies exceeding 100 TeV. Due 
to their very small flux at these energies, cosmic 
rays can only be detected indirectly via the ex- 
tensive air showers (EAS) they induce in the at- 
mosphere. Hence, in order to infer the mass and 
energy of the primary cosmic rays as well as to 
gain insight into the physics of the high-energy 
collisions, one has to rely on detailed simulations 
of air showers. For this purpose, hadronic mul- 
tiparticle production has to be simulated at en- 
ergies exceeding by far those accessible at man- 
made accelerators and in phase space regions not 
covered in collider experiments. The secondary 
muons produced in these collisions which can be 
detected directly on ground, are a tracer of high- 
energy hadronic interactions. For this reason, the 
number of muons produced in simulated air show- 
ers depends strongly on the adopted hadronic in- 
teraction models [TJ. In addition, muons are a 
sensitive indicator of the primary composition, 
as, for example, iron showers produce about 40% 
more muons than protons. 

In this work we will employ air shower univer- 
sality in order to determine the number of muons 
in air showers detected with Auger, using sev- 
eral independent methods. Air shower univer- 
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sality |2|3j states that the average properties of 
air showers can be completely characterized by 
the primary energy, stage of shower evolution, 
and overall muon normalization, when viewing 
the shower around or after shower maximum. In 
particular, the electromagnetic particle content of 
the shower is already determined by its energy 
and stage of evolution. Physically, this is due 
to the fact that the electromagnetic cascade in- 
volves dozens of particle generations and of order 
10 9 — 10 10 particles, which washes out the details 
of high-energy hadronic interactions. 

In order to infer the muon density of showers 
at a given energy, one has to rely on an energy 
scale of the surface detector. Using the univer- 
sality of the electromagnetic contribution to the 
signal measured by the Auger surface detectors, 
5(1000), we are able to determine an energy scale 
independent of the fluorescence energy measure- 
ment. The muon density at ground which we infer 
using the surface detector and, independently, us- 
ing hybrid events, does not rely on assumptions 
on the primary cosmic ray composition. This al- 
lows for a direct test of the predictions of hadronic 
interaction models. 

2. Parameterization of the average surface 
detector signal using universality 

Universality features of the longitudinal profile 
of showers have been studied by several authors. 
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Here we exploit shower universality features to 
predict the surface detector signal expected for 
Auger Cherenkov tanks from the electromagnetic 
and muonic shower components at 1000 m from 
the shower core. In the following, we give a brief 
outline of the parameterization. For a detailed 
description, see [3]. 

A library of proton and iron showers cover- 
ing the energy range from 10 17 to 10 20 eV and 
zenith angles between 0° and 70° was generated 
with CORSIKA 6.5 [4] and the hadronic inter- 
action models QGSJET 11.03 [5 and FLUKA 
[6]. For comparison, a smaller set of showers 
was simulated with the combinations QGSJET 
II.03/GHEISHA [7] and SIBYLL 2.1/FLUKA 
8 9 J . Seasonal models of the Malargiie molecular 
atmosphere were used [TO] . The detector response 
is calculated using look-up tables derived from a 
detailed GEANT4 simulation of the Auger sur- 
face detectors, which has been shown to match 
the data well [TTj . Note that since the Auger 
detectors are calibrated against the cosmic ray 
muon background, we only rely on the correct 
simulation of the signal generated by electromag- 
netic particles relative to that of the muons. 

Air shower universality states that the electro- 
magnetic component of the predicted surface de- 
tector signal ,5(1000) at the lateral distance of 
1000 m depends only on primary energy and the 
stage of shower evolution, and not on the pri- 
mary particle and hadronic model assumed. We 
define the electromagnetic component of the sig- 
nal as that of electrons, positrons, and gamma 
rays excluding muon decay products, and mea- 
sure shower evolution in terms of the distance 
between the shower maximum and the ground 
DX = A groun d — A max , measured along the 
shower axis. In fact, universality is slightly vi- 
olated, and the electromagnetic contribution to 
S'(IOOO) of proton and iron showers differ by 
about 12% [3], see Fig. [TJ The dependence on 
the hadronic model is much smaller, about 5%. 

The muon contribution to 5(1000), which in- 
cludes muon decay products in our definition, 
depends strongly on the primary particle (40% 
difference between proton and iron) and the 
hadronic model (Fig.[T]). However, the functional 
dependence on DX is universal. 
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Figure 1. Upper panel: Simulated electromagnetic 
shower plane signals (excluding muon decay prod- 
ucts) at 1000 m for proton (red dots) and iron showers 
(blue circles) at 10 19 eV [6 = 0-60°; using QGSJetll 
and Fluka) as a function of DX. Lower panel: Simu- 
lated muon signals at 1000 m (including muon decay 
products) vs. DX for the same simulated showers, 
and proton and iron showers simulated using Sibyll 
(crosses). The lines show joint Gaisser-Hillas type 
parameterizations of all primaries and models. 



After accounting for geometrical effects such 
as the projected tank surface area, the electro- 
magnetic shower signal (averaged over proton and 
iron) is parameterized as function of the energy 
E, distance to shower maximum DX, and zenith 
angle 9. The difference between proton and iron 
shower profiles is included in the calculation of the 
systematic uncertainties later on. Similarly, the 
universal shape of the muon signal profile is pa- 
rameterized simultaneously using a Gaisser-Hillas 
function for all models and primaries (Fig. [I}, 
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leaving only normalization factors free. The ex- 
pected detector signal at 1000 m can then be writ- 
ten as 

S'param(-E') 9, X max , A^J = S cm (E, 9, DX(X mSLX )) 

+N„ S^ GSII ^(10EeV,9,DX(X max )), (1) 

where ,p is the parametrized muon signal 

adopting the normalization for proton-QGSJET 
II at 10 EeV, so that 7V M is the number of muons 
relative to that of QGSJET-II proton showers at 
10 EeV. 

3. Measuring from the surface detector 

Within the current statistics, the arrival direc- 
tions of high-energy cosmic rays around 10 EeV 
are found to be isotropic, allowing us to apply 
the constant intensity method to determine the 
muon signal contribution. Dividing the surface 
detector data into equal exposure bins, the num- 
ber of showers with 5(1000) greater than a given 
threshold S param (lOEeV,0,(X max ),2V M ), should 
be the same for each bin. 



dN n 



dsin 



const. (2) 
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Using the independently measured mean depth of 
shower maximum (V max ) [13] . the only remaining 
free parameter in Eq. @ is the relative number of 
muons Np. For a given energy E, is adjusted 
to obtain a flat distribution of events in sin 2 9. 

More precisely, scanning through a range of 
values, we calculate the x 2 /dof of the event his- 
togram relative to a flat distribution in sin 2 9. We 
then fit a two-branch parabola around the mini- 
mum of x 2 (V M ), which results in a best-fit N^m 
and its asymmetric error, a N . 

The sensitivity of this method to the muon 
number parameter in Eq. JT]) is illustrated in 
Fig. O Clearly, = 1 (top histogram in the 
lower panel of Fig. \2§ is ruled out by the data. 
The best description of the data above 10 EeV 
requires N^fn = 1.60. However, we have to take 
into account shower-to-shower fluctuations and 
the finite resolution of 5(1000) which lead to a 
slightly biased iV M measurement in this method. 

In order to estimate this bias, we have sim- 
ulated a large number of Monte Carlo realiza- 
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Figure 2. Upper panel: the signal parameterization 
Eq. ([TJ vs. sec for different iV^ (black/solid— 1.6, 
red/dashed— 1.0, blue/dotted— 2.0). Lower panel: 
histograms of number of Auger events above the pa- 
rameterized signal in equal exposure bins, obtained 
for the same iV^ as shown in the upper panel. The 
black histogram is for iV M = 1.6, the best-fit value 
found in the data (see text). 

tions of SD data sets of the same size as the 
Auger data set, distributed according to the ob- 
served cosmic ray spectrum and assuming dif- 
ferent primary compositions (pure proton, iron, 
or mixed composition), and different "true" N^. 
For each shower, X max is obtained from distri- 
butions which closely match the observed X max 
distributions in the fluorescence data. The elec- 
tromagnetic and muon signals are fluctuated ac- 
cording to the model predictions [3]. Note that 
the magnitude of fluctuations in V max and 
are robust predictions which only depend on the 
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primary particle. Eq. JT]) is then used to calcu- 
late the signal at 1000 m from the shower core, 
5(1000), which is smeared with an experimental 
reconstruction accuracy [12] . A realistic trigger 
probability is also applied. 

For each simulated data set, we then apply ex- 
actly the same analysis as used on the real data. 
We found that is systematically overestimated 
by 3—6%. This bias is due to a combined action 
of shower-to-shower fluctuations and detector res- 
olution which depends on the signal size. We 
have found that the bias in is independent 
of the "true" value of iV^ adopted, and primar- 
ily depends on the size of reconstruction uncer- 
tainties on 5(1000) and the primary composition. 
Fig. [3] shows the relative bias on as a function 
of primary energy, for different assumed compo- 
sitions and increased/decreased detector resolu- 
tion. Clearly, the bias is quite robust to rather 
extreme changes in those assumptions. 

We then subtract the bias expected for a mixed 
composition (black dots in Fig. [3]) from the mea- 
sured Np fit, obtaining a corrected of: 

JV M (10 EeV) = 1.53±g;g? (stat.)±°;^ (syst.) (3) 

Three main sources of systematics enter in the 
determination of using the constant intensity 
method: the uncertainty on the electromagnetic 
signal due to universality violation (Sec. [2]); the 
uncertainty on (X max ) at this energy; and the 
uncertainty on the bias of determined from 
the Monte Carlo simulations. 

In order to quantify the first two systematics, 
we measure 7V M using an electromagnetic signal 
rescaled by ±6%, bracketing the observed univer- 
sality violation, and measuring with varying 
(I mffl ), adopting the systematic and statistical 
uncertainty reported by Auger |13) . For the un- 
certainty on the bias, we adopt ±3%. Summing 
the different contributions in quadrature, we ob- 
tain the total systematic error on quoted in 
Eq. ©. 

Once N,j, at 10 EeV is measured, and using the 
measured mean depth of shower maximum, the 
signal size at 9 = 38° can be calculated: 

538(10 EeV) = 38.91^ (stat.)tlj (sys.) VEM (4) 
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Figure 3. Relative bias on iV M measured with the 
constant intensity method, determined from Monte 
Carlo simulations of Auger data sets. The differ- 
ent points show various assumptions on the primary 
composition and the experimental resolution of the 
5(1000) reconstruction. 



Within the uncertainties, this value of S38 is con- 
sistent with the energy scale from fluorescence 
detector measurements, whose systematic uncer- 
tainty (22% [14]) is dominated by the uncertainty 
on the fluorescence yield. Eq. ^ corresponds to 
assigning showers a ~ 26% higher energy than 
done in the fluorescence detector-based Auger 
shower reconstruction (E = 1.26 2?fd)- 

4. Measuring from hybrid events 

The Pierre Auger Observatory has the unique 
capability of measuring hybrid events which have 
been simultaneously detected by the surface de- 
tector array and fluorescence telescopes. For each 
of these events, a calorimetric measurement of the 
energy E and depth of shower maximum X max 
are available, in addition to the signal 5(1000) 
from the surface array. Hence, we can use our 
parameterization of the universal electromagnetic 
signal S em (E, A max , 9) to determine the electro- 
magnetic contribution to 5(1000). The remainder 
of the signal is attributed to muons, and we can 
determine a muon normalization for each event: 



5(1000) - S cm (E, X 

max 1 u J 

Sn(0. ^raax 



(5) 
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Figure 4. Reconstructed muon tank signal contri- 
bution as a function of shower evolution (distance to 
ground DX) for hybrid events with 9 < 60° (black 
points) and 60° < 9 < 68° (grey points), using an 
energy scale of E = 1.26-Efd- The muon signal from 
simulated proton (red) and iron (blue) showers (using 
QGSJetll) is also shown. 



where is again the reference muon signal 
(proton-QGSJetll at 10 EeV). 

We select high-quality hybrid events for which 
the shower maximum is in the field of view of a 
telescope, 8 < 60°, and the Mie scattering length 
is measured. Furthermore, the Cherenkov light 
fraction is limited to less than 50%, and we apply 
all cuts used in the measurement of the elongation 
rate [13] which have been carefully chosen in order 
to ensure an unbiased V max distribution of the air 
showers. The surface detector event has to satisfy 
the T5 selection cuts which are also applied in 

The average muon normalization in hybrid 
events at 10 EeV is found to be 

N »\e=i.2 6 e fd = 1-49 ±0.05, and 

N^\ e=Efd = 1.84 ±0.07, (6) 

in excellent agreement with the constant inten- 
sity method results, when compared at the same 
energy scale. 

A similar study has been performed for inclined 
hybrid events (60° < 6 < 68°). While the statis- 
tics are limited, no subtraction of S om is necessary 
for these events, as the electromagnetic part has 



Figure 5. Comparison of the results on at 10 EeV 
from different methods. The black dot shows the con- 
stant intensity method result with statistical (diag- 
onal line) and systematic errors (vertical/horizontal 
lines). The blue line with shaded band, and red lines 
indicate the result from vertical and inclined hybrid 
events, respectively, both with statistical errors. The 
vertical dashed lines indicate the systematic error on 
the fluorescence energy scale [14] . 



been attenuated in the large depth of atmosphere 
traversed. Thus, inclined hybrid events yield a 
clean measurement of N^, and good agreement 
between muon numbers of the inclined and the 
vertical data sets is found. 

Fig. 2] shows the muon signal derived from the 
hybrid events as a function of shower evolution 
(DX). While the muon signal is clearly higher 
than that predicted in the simulations, the be- 
havior as a function of DX follows the prediction 
very well, which serves as a consistency test of 
our method. 

In Fig. [5] we compare the results of the different 
methods applied for inferring the muon density at 
1000 m from the shower core. The relative num- 
ber of muons is shown as function of the adopted 
energy scale with respect to the Auger fluores- 
cence detector energy reconstruction. Only the 
constant intensity method yields an independent 
measurement of the energy scale. Good agree- 
ment between the presented methods is found, 
when compared at E ~ 1.26 Epjj. 
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5. Discussion 

Assuming air shower universality, which was 
tested with simulations of various primaries and 
hadronic interaction models, we have determined 
the muon density at 10 EeV and the energy scale 
with which the data of the Pierre Auger Observa- 
tory can be described self-consistently. The num- 
ber of muons measured in data is about 1.5 (1.7) 
times larger than that predicted by QGSJET II 
(Sibyll) for proton showers, with statistical and 
systematic errors of the order of 0.1. This value 
of Nu is even slightly larger than that predicted 
for iron showers, by a factor of 1.1 and 1.2 for 
QGSJetll and Sibyll, respectively. A similar ex- 
cess has been found at energies from 3 EeV to 
50 EeV. 

Consistent results were obtained with several 
analysis methods on independent data sets: the 
constant intensity method (Sec. [3]) uses only the 
surface detector data and is independent of the 
fluorescence energy scale. Vertical hybrid events 
(Sec. S]) allow for a direct and unbiased measure- 
ment of for each event, albeit depending on 
the flourescence energy. The signal in inclined hy- 
brid events (9 > 60°) is purely muonic, and the 
resulting independent of the subtraction of 
the electromagnetic signal. The fact that all three 
methods agree within errors for an energy scale 
of 1.26-Efd indicates that increasing the number 
of muons by 50% (relative to proton-QGSJetll) 
yields a consistent description of the Auger data. 

If interpreted in terms of cosmic ray compo- 
sition, the derived muon density would corre- 
spond to nuclei heavier than iron, which is clearly 
at variance with the measured V max values at 
this energy. The discrepancy between air shower 
data and simulations reported here is qualita- 
tively similar to the inconsistencies found in com- 
position analyses of previous detectors [15116117] . 
This points towards a deficiency of hadronic mod- 
els in predicting the number of muons produced 
in air showers at large distances from the shower 
core. However, it is important to note that the 
method presented here relies on a correct descrip- 
tion of the electromagnetic shower component in 
the simulations. 
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